GEO UNIV’R Tunisie 2024

Faire des cartes thématiques avec R

Author

Nicolas Lambert, Elina Marveaux, Ronan Ysebaert

Published

April 22, 2024

Note

Ce support est très largement inspiré du manuel Cartographie avec R réalisé par Timothée Giraud et Hugues Pécout. Et on les en remercie chalereusement 🙏

Pakages utilisés dans cette session.

  • sf
  • mapsf
  • cartogram (facultatif)
  • leaflet (facultatif)
  • mapview (facultatif)
  • wbstats (facultatif)

Initiez un nouveau projet

  • Créer un nouveau repertoire de travail.
  • Ouvrez Rstudio.
  • Créez un nouveau projet et placez-le dans votre dossier.
  • Créez un nouveau script R (ou un document Quarto)

Si vous n’avez pas ces packages listés plus haut, vous pouvez les installer en tapant la ligne suivante dans la console.

install.packages(c('sf', 'mapsf', 'cartogram', 'leaflet', 'mapview', 'wbstats'))

Pour gagner du temps, seuls sf et mapsf sont absolument nécéssaires. Vous pouvez donc installer uniquement ces deux là.

install.packages(c('sf', 'mapsf'))

Dans cette séquence, nous travaillons à l’échelle des pays africains. Les données statistiques sont issues de la base de donnée de la banque mondiale.

Les données sont disponibles au téléchargement ici.

Téléchargez-les, dézippez-les et placez-les dans un répertoire data.

Tout est prêt 😎

Import et mise en forme des données

1 - Import des géométries

library("sf")

On dispose d’un fichier geipackage contenant plusieurs géométries

st_layers("data/afrique.gpkg")
Driver: GPKG 
Available layers:
               layer_name     geometry_type features fields crs_name
1 ne_10m_populated_places             Point     1305    137   WGS 84
2            ne_10m_roads Multi Line String     6761     31   WGS 84
3        ne_10m_railroads Multi Line String      758     12   WGS 84
4                  africa           Polygon       55      4   WGS 84
5                   world           Polygon      209      4   WGS 84

On importe la couche correspondant aux pays africains et on la reprojette en projection Pseudo-Mercator.

geom <- st_transform(st_read("data/afrique.gpkg", layer = "africa"),"epsg:3857")
Simple feature collection with 4 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2702989 ymin: -2329942 xmax: 6434678 ymax: 1782119
Projected CRS: WGS 84 / Pseudo-Mercator
  ISO3            NAMEen             NAMEfr region
1  GNQ Equatorial Guinea Guinée-Équatoriale Africa
2  COM           Comoros            Comores Africa
3  CPV        Cape Verde           Cap-Vert Africa
4  MUS         Mauritius            Maurice Africa
                            geom
1 MULTIPOLYGON (((1260368 241...
2 MULTIPOLYGON (((4883099 -13...
3 MULTIPOLYGON (((-2660139 17...
4 MULTIPOLYGON (((6424679 -22...

2 - Import des données statistiques issues de la banque mondiale.

data <- read.csv("data/worldbank_africa.csv")
          id  type                                      description
1     region quali                                        subregion
2        pop stock                                Population, total
3        urb stock                                 Urban population
4      rural stock                                 Rural population
5   ruralpct ratio         Rural population (% of total population)
6        gdp stock                                GDP (current US$)
7      gdppc ratio                     GDP per capita (current US$)
8        esp ratio          Life expectancy at birth, total (years)
9       elec ratio          Access to electricity (% of population)
10    forest stock                             Forest area (sq. km)
11 forestpct ratio                     Forest area (% of land area)
12  articles stock        Scientific and technical journal articles
13  internet ratio Individuals using the Internet (% of population)
14   servers ratio   Secure Internet servers (per 1 million people)

3 - Jointure

africa <-  merge(
  x = geom[,"ISO3"],  
  y = data,  
  by.x = "ISO3",
  by.y = "id",
  all.x = TRUE   
)

Et quelques couches additionnelles

world <- st_transform(st_read("data/afrique.gpkg", layer = "world"),"epsg:3857")
places <- st_transform(st_read("data/afrique.gpkg", layer = "ne_10m_populated_places"),"epsg:3857")
rail <- st_transform(st_read("data/afrique.gpkg", layer = "ne_10m_railroads"),"epsg:3857")
roads <- st_transform(st_read("data/afrique.gpkg", layer = "ne_10m_roads"),"epsg:3857")

1 - cartographie avec sf

Grâce au package sf, la fonction plot() appliquée à un Spatial*DataFrame permet de l’afficher comme n’importe quel graphique.

Comportement par defaut

Par defaut, la fonction renvoie une image pouvant contenir jusqu’à 10 cartes, chaque carte correspondant à une colonne.

plot(africa)

Si il y a plus de 10 colonnes dans le jeu de données, alors on ajoute la paramètre max.plot = ncol(africa) - 1 (le nombre de colonnes total moins la colonne de géometries).

plot(africa, max.plot = ncol(africa) - 1)

On perçoit que les couleurs ne sont pas choisies totalement au hasard. Nous verrons un peu plus tard comme cela fonctionne.

Si on souhaite afficher simplement les géométries, on utilise la fonction st_geometry()

plot(st_geometry(africa))

Ou alors, en sélectionnant le champ contenant les géométries.

plot(africa$geometry)

Marges en emprise

Dans les affichages précédents, on constate qu’il y a beaucoup d’espace perdu autour de la carte. On peut modifier les marges en modifiant les parametres graphiques. mar permet de mettre les marges à zero. xaxs='i', yaxs='i' permettent de supprimer les espaces residuels.

par(mar = c(0, 0, 0, 0),  xaxs='i', yaxs='i', bg = "#F1F3F5") # c(bottom, left, top, right))
plot(st_geometry(africa))

Quand on travaille dans un document quarto, il est également possible de modifier la taille des figures avec les parametres fig-height et fig-width. Soit dans le yaml en haut du document. Soit directement dans le chunk.

Ici, on détermine le bon ratio entre la hauteur et la largeur.

bb <- st_bbox(africa)
ratio <- (bb$xmax - bb$xmin) / (bb$ymax - bb$ymin) 
height <- 6
width <- height * ratio
width
    xmax 
6.355375

Puis on écrit

#| fig-heigth: 6
#| fig-width: 6.36

Et on obtient une carte qui prend bien tout l’espace souhaité.

par(mar = c(0, 0, 0, 0), bg = "#F1F3F5", xaxs='i', yaxs='i')
plot(st_geometry(africa))

Enfin, notez qu’avec xlim et ylim, vous pouvez cadrer la carte sur une emprise particulière. Par exemple ici, on affiche le fond de carte du monde mais uniquement avec l’emprise de l’Afrique.

par(mar = c(0, 0, 0, 0), xaxs='i', yaxs='i')
bb <- st_bbox(africa)
plot(st_geometry(world), xlim = c(bb$xmin, bb$xmax), ylim = c(bb$ymin, bb$ymax))

Paramétrer le style

De nombreux paramètres permettent de personnaliser le style de la carte. Nous en citons ici quelques exemples.

  • Couleur de fond : col
  • Couleur de contour : border
  • Epaisseur des tracés : lwd
par(mar = c(0, 0, 0, 0), bg = "#F1F3F5", xaxs='i', yaxs='i')
plot(st_geometry(africa), col = "#5b89a3", border = "white", lwd = 0.5)

NB : Il y a 657 noms de couleurs disponibles dans R. Pour les afficher, vous pouvez taper colors()

cols <- colors()
head(cols, 20)
 [1] "white"         "aliceblue"     "antiquewhite"  "antiquewhite1"
 [5] "antiquewhite2" "antiquewhite3" "antiquewhite4" "aquamarine"   
 [9] "aquamarine1"   "aquamarine2"   "aquamarine3"   "aquamarine4"  
[13] "azure"         "azure1"        "azure2"        "azure3"       
[17] "azure4"        "beige"         "bisque"        "bisque1"

Avec lty, on peut également changer le type de traits

par(mar = c(0, 0, 0, 0), bg = "#F1F3F5", xaxs='i', yaxs='i')
plot(st_geometry(roads), lty = 3)

Avec pch, on peut choisir le s de symbole. Avec cex, on determine sa taille.

par(mar = c(0, 0, 0, 0), bg = "#F1F3F5", xaxs='i', yaxs='i')
plot(st_geometry(places), pch = 17, col = "red", cex = 1)

Supperposer des couches

Pour supproposer des couches, vous pouvez simplement utiliser le paramètre add = TRUE (sous reserve que les couches soient dans le même système de coordonnées)

par(mar = c(0, 0, 0, 0), bg = "#F1F3F5")
plot(st_geometry(africa), col = "#5b89a3", border = "white", lwd = 0.5)
plot(st_geometry(roads), col = "red", add = TRUE)
plot(st_geometry(places), pch = 19, col = 'black', cex = 0.5, add = TRUE)

Autres éléments

Il est possible d’ajouter d’autres éléments à la mise en page.

  • Les axes : axes = TRUE/FALSE
  • Graticule : graticule = TRUE
  • Titre = main = "Hello"

Par exemple :

par(mar = c(0, 0, 2, 0))
plot(st_geometry(africa), col= "white", axes = TRUE, graticule = TRUE, main = "Hello")

Effet d’ombrage

Avec R, il est aisé de translater une géoémétrie pour créer un effet d’ombrage.

par(mar = c(0, 0, 0, 0))
plot(st_geometry(africa) + c(50000,-50000), col = "#827e6c80", border = NA)
plot(st_geometry(africa) , col = "#5B89A3", border = NA, add = TRUE)

Astuce

On peut faire varier la transparence d’une couleur au fromat hexadecimal en ajouter un nombre de 00 à 99 à la fin du code. Par exemple #827e6c60 applique une opacité de 60% à la couleur #827e6c.

Cartographie thématique

Avec le package sf, il est (un peu) possible de réaliser des cartes thématiques. Rappelez-vous les couleurs de tout à l’heure. Dans le cas où un seul attribut est sélectionné, une légende est attribuée par défaut à côté de la carte. Ici, une donnée qualitative.

plot(africa["region"])

Le positionnement de la légende peut être défini par le paramètre key.pos (1 = dessous, 2 = gauche, 3 = dessus and 4 = droite). Sa taille peut également être modifiée avec les paramètres key.width et key.length

plot(africa["region"], key.pos = 1, key.length = 1)

Si on souhaite cartographier une variable quantitative, la palette par defaut est différente.

plot(africa["pop"])

Grace au paramètre breaks, il est possible de donner ses propres classes de valeur ou de donner une méthode de discrétisation (méthodes du package classInt).

plot(africa["pop"], breaks = "jenks")

Exports

Ici, on a affiché toutes les cartes dans le document. Mais on peut également choisir de les construire dans un format donné (pdf, svg, png, ps, etc.), ce qui peut être utile pour les retravailler dans un logiciel de DAO. Par exemple, on peut écrire :

svg("my_plot.svg")
plot(st_geometry(africa))
dev.off()

Bilan et limites

Que retenir ?

Les fonctions de cartographies dans les fonctionnalités de base de sf sont très limitées. On ne peut pas, par exemple, dessiner des symboles proportionnels et leur légende associée. Si on veut aller plus loin, on a besoin d’un package spécialisé en représentations cartographiques. C’est à ce besoin que répond le package mapsf.

2 - le package mapsf

mapsf permet de créer la plupart des types de carte utilisés habituellement en cartographie statistique (cartes choroplèthes, typologies, symboles proportionnels ou gradués…).

Pour chaque type de carte, plusieurs paramètres permettent de personnaliser la représentation cartographique. Ces paramètres sont les mêmes que ceux que l’on retrouve dans les logiciels de SIG ou de cartographie usuels. mapsf est le successeur du package cartography. Si vous avez l’habitude d’utiliser ce dernier, nous vous conseillons vivement d’utiliser dorénavant mapsf.

On charge le package

library("mapsf")

Documentation et supports

De nombreux documents permettent de prendre en main ce package.

Et surtout, il faut aller voir dans la documentation du package directement dans RStudio. Vous y decouvrirez une magnifique vignette.

Afficher un fond de carte

La fonction mf_map() est la fonction centrale du package mapsf. Elle permet de réaliser la plupart des représentations usuelles en cartographie. Par défaut, elle permet d’afficher un fond de carte.

mf_map(x = africa, border = NA, col = NA)
mf_map(x = world, border = "white", col = "#CCCCCC50", lwd = 0.5, add = TRUE)
mf_map(x = africa, border = "white", col = "#6893d9", lwd = 0.5, add = TRUE)
mf_map(x = places, pch = 20, cex = .7, col = "darkred", add = TRUE)
mf_title(txt = "L'Afrique")

Symboles proportionnels

Pour représenter une donnée quantitative absolue (i.e. donnée de stock), on utilise la fonction mf_map avec le paramètre type = "prop"

mf_map(x = africa, border = "white", lwd = 0.5)
mf_map(x = africa,
       var = "pop",
       type = "prop",
       border = "white",
       col = "#FF000080",
       leg_title = "Nombre d'habitants\nen 2020",
       inches   = 0.4)
mf_title(txt = "Population totale")

Carte choroplèthe

Carte de typologie

Combinaisons

Habillage et mise en page

A vous de jouer

Réaliser une carte sur un indicateur de votre choix. Si besoin, vous pouvez aller chercher d’autres indicateurs grace au package wbstats

library(wbstats)

Par exemple

wb_search(pattern = "Bird")
# A tibble: 2 × 3
  indicator_id   indicator                  indicator_desc                      
  <chr>          <chr>                      <chr>                               
1 EN.ANM.THRD.NO Animal species, threatened Animal species are mammals (excludi…
2 EN.BIR.THRD.NO Bird species, threatened   Birds are listed for countries incl…
oiseaux_menaces <- wb_data("EN.BIR.THRD.NO", start_date = 2018, end_date = 2018)
head(oiseaux_menaces)
# A tibble: 6 × 9
  iso2c iso3c country              date EN.BIR.THRD.NO unit  obs_status footnote
  <chr> <chr> <chr>               <dbl>          <dbl> <chr> <chr>      <chr>   
1 AW    ABW   Aruba                2018              2 <NA>  <NA>       <NA>    
2 AF    AFG   Afghanistan          2018             16 <NA>  <NA>       <NA>    
3 AO    AGO   Angola               2018             32 <NA>  <NA>       <NA>    
4 AL    ALB   Albania              2018              8 <NA>  <NA>       <NA>    
5 AD    AND   Andorra              2018              3 <NA>  <NA>       <NA>    
6 AE    ARE   United Arab Emirat…  2018             13 <NA>  <NA>       <NA>    
# ℹ 1 more variable: last_updated <date>
africa_birds <-  merge(
  x = geom[,"ISO3"],  
  y = oiseaux_menaces,  
  by.x = "ISO3",
  by.y = "iso3c",
  all.x = TRUE   
)
mf_map(x = africa, border = NA, col = NA)
mf_map(x = world, border = "white", col = "#CCCCCC50", lwd = 0.5, add = TRUE)
mf_map(x = africa, border = "white", col = "#6893d9", lwd = 0.5, add = TRUE)
mf_map(x = africa_birds,
       var = "EN.BIR.THRD.NO",
       type = "prop",
       symbol = "square",
       border = "white",
       col = "#FF000080",
       leg_title = "Nombre d'oiseaux\nmenacés en 2018",
       inches   = 0.3)

3 - Cartographie avancée

Cartograms

Il existe plusieurs méthodes pour réaliser des cartogrammes.

1. Cartogrammes de Dorling

library(cartogram)
pop2020_dorling <- cartogram_dorling(
  africa[!is.na(africa$pop),],
  weight = "pop",
  k=2.5
  )
mf_map(africa, col = "white", border= NA)
mf_map(pop2020_dorling, col = "#5B89A3", border= "white", add = TRUE)
mf_label(
  x = pop2020_dorling[order(pop2020_dorling$pop, decreasing = TRUE), ][1:10,],
  var = "name",
  col = "#5B89A3",
  overlap = FALSE, lines = FALSE,
  halo = TRUE,
  r = .15
)
mf_title("Population totale - Cartogramme de Dorling")

2. Les cartogrammes non continus

afr_ncont <- cartogram_ncont(x = africa, weight = "pop", k = 1.2)
mf_map(africa, border = "white", lwd = 0.5,)
mf_map(afr_ncont, col = "#5B89A3", border= "white", add = TRUE)
mf_title("Population en Afrique - Cartogramme de Olson")

3. Cartogrammes continus

africa[is.na(africa$pop),"pop"] <- 1
afr_cont <- cartogram_cont(x = africa,
                           weight = "pop",
                           itermax = 30)
mf_map(afr_cont, col = "#5B89A3", border= "white", add = FALSE)
mf_title("Population en Afrique - Cartogramme de Dougenik")
mf_inset_on(africa, cex = .2, pos = "topleft")
mf_map(africa, lwd = .5, border = "white")
mf_inset_off()

Carroyage

Discontinuités

???

Lissages

???

mapview

library(mapview)
mapview(africa) + mapview(st_centroid(africa))
Warning: st_centroid assumes attributes are constant over geometries

Leaflet

Leaflet est un package basé sur le JavaScript, permettant de faire de la cartographie interactive.

library(leaflet)

Réalisation d’une primière carte simple

m = leaflet() %>% addTiles()
m

Zoom sur une localisation précise.

sfax <- c(10.760034694759957, 34.7407779744004)
m2 <- leaflet() %>% setView(lng = sfax[1], lat = sfax[2], zoom = 12) %>% 
  addTiles() 
m2

Ajout de géométries

africa_wgs84 <- st_transform(africa, 4326)
popup <- paste0("<b>",africa_wgs84$name,"</b><br/><b>Population: </b>", 
                africa_wgs84$pop)
m3 = leaflet() %>% 
  addTiles() %>% 
  addPolygons(data=africa_wgs84, weight = 2, fillColor = "yellow", popup= popup) %>%         
  addMarkers(data = st_centroid(africa_wgs84)) %>%  addMiniMap(position = "bottomright")
Warning: st_centroid assumes attributes are constant over geometries
m3

Au delà de R

Avec Quarto, il est également possible de créer dans Rstudio, des cartes thématiques interactives pour le

Pour celles et ceux qui le souhaitent, c’est ce que nous ferons vendredi matin.